{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center>\n",
    "<img src=\"../../img/ods_stickers.jpg\">\n",
    "## Открытый курс по машинному обучению\n",
    "<center>Авторы материала: Павел Нестеров, Юрий Кашницкий\n",
    "    \n",
    "Материал распространяется на условиях лицензии [Creative Commons CC BY-NC-SA 4.0](https://creativecommons.org/licenses/by-nc-sa/4.0/). Можно использовать в любых целях (редактировать, поправлять и брать за основу), кроме коммерческих, но с обязательным упоминанием автора материала."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>Тема 4. Линейные модели классификации и регрессии\n",
    "## <center>Часть 1. Линейная регрессия"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Метод наименьших квадратов\n",
    "\n",
    "Рассказ про линейные модели мы начнем с линейной регрессии. В первую очередь, необходимо задать модель зависимости объясняемой переменной $y$ от объясняющих ее факторов, функция зависимости будет линейной: $y = w_0 + \\sum_{i=1}^m w_i x_i$. Если мы добавим фиктивную размерность $x_0 = 1$ для каждого наблюдения, тогда линейную форму можно переписать чуть более компактно, записав свободный член $w_0$ под сумму: $y = \\sum_{i=0}^m w_i x_i = \\textbf{w}^{\\text{T}} \\textbf{x}$. Если рассматривать матрицу наблюдения-признаки, у которой в строках находятся примеры из набора данных, то нам необходимо добавить единичную колонку слева. Зададим модель следующим образом:\n",
    "\n",
    "$$\\large \\textbf y = \\textbf{X} \\textbf w + \\epsilon,$$\n",
    "\n",
    "где\n",
    "- $\\textbf{y} \\in \\mathbb{R}^n$ – объясняемая (или целевая) переменная;\n",
    "- $\\textbf{w}\\in \\mathbb{R}^{m+1}$ – вектор параметров модели (в машинном обучении эти параметры часто называют весами);\n",
    "- $\\textbf{X}$ – матрица наблюдений и признаков размерности $n$ строк на $m + 1$ столбцов (включая фиктивную единичную колонку слева) с полным рангом по столбцам: $\\text{rank}\\left(\\textbf{X}\\right) = m$;\n",
    "- $\\epsilon$ – случайная переменная, соответствующая случайной, непрогнозируемой ошибке модели.\n",
    "\n",
    "Можем выписать выражение для каждого конкретного наблюдения\n",
    "\n",
    "$$\\large \n",
    "y_i = \\sum_{j=0}^m w_j X_{ij} + \\epsilon_i$$\n",
    "\n",
    "Также на модель накладываются следующие ограничения (иначе это будет какая то другая регрессия, но точно не линейная):\n",
    "- матожидание случайных ошибок равно нулю: $\\forall i: \\mathbb{E}\\left[\\epsilon_i\\right] = 0$;\n",
    "- дисперсия случайных ошибок одинакова и конечна, это свойство называется <a href=\"https://ru.wikipedia.org/wiki/%D0%93%D0%BE%D0%BC%D0%BE%D1%81%D0%BA%D0%B5%D0%B4%D0%B0%D1%81%D1%82%D0%B8%D1%87%D0%BD%D0%BE%D1%81%D1%82%D1%8C\">гомоскедастичностью</a>: $\\forall i: \\text{Var}\\left(\\epsilon_i\\right) = \\sigma^2 < \\infty$;\n",
    "- случайные ошибки не скоррелированы: $\\forall i \\neq j: \\text{Cov}\\left(\\epsilon_i, \\epsilon_j\\right) = 0$.\n",
    "\n",
    "Оценка $\\widehat{w}_i$ весов $w_i$ называется *линейной*, если\n",
    "\n",
    "$$\\large \\widehat{w}_i = \\omega_{1i}y_1 + \\omega_{2i}y_2 + \\cdots + \\omega_{1n}y_n,$$\n",
    "\n",
    "где $\\forall\\ k\\ \\omega_{ki}$ зависит только от наблюдаемых данных $\\textbf{X}$ и почти наверняка нелинейно. Так как решением задачи поиска оптимальных весов будет именно линейная оценка, то и модель называется *линейной регрессией*. Введем еще одно определение. Оценка $\\widehat{w}_i$ называется *несмещенной* тогда, когда матожидание оценки равно реальному, но неизвестному значению оцениваемого параметра:\n",
    "\n",
    "$$\\large \\mathbb{E}\\left[\\widehat{w}_i\\right] = w_i$$\n",
    "\n",
    "Один из способов вычислить значения параметров модели является <b>метод наименьших квадратов</b> (МНК), который минимизирует среднеквадратичную ошибку между реальным значением зависимой переменной и прогнозом, выданным моделью:\n",
    "\n",
    "$$\\large \\begin{array}{rcl}\\mathcal{L}\\left(\\textbf{X}, \\textbf{y}, \\textbf{w} \\right) &=& \\frac{1}{2n} \\sum_{i=1}^n \\left(y_i - \\textbf{w}^{\\text{T}} \\textbf{x}_i\\right)^2 \\\\\n",
    "&=& \\frac{1}{2n} \\left\\| \\textbf{y} - \\textbf{X} \\textbf{w} \\right\\|_2^2 \\\\\n",
    "&=& \\frac{1}{2n} \\left(\\textbf{y} - \\textbf{X} \\textbf{w}\\right)^{\\text{T}} \\left(\\textbf{y} - \\textbf{X} \\textbf{w}\\right)\n",
    "\\end{array}$$\n",
    "\n",
    "Для решения данной оптимизационной задачи необходимо вычислить производные по параметрам модели, приравнять их к нулю и решить полученные уравнения относительно $\\textbf w$ (матричное дифференцирование неподготовленному читателю может показаться затруднительным, попробуйте расписать все через суммы, чтобы убедиться в ответе):\n",
    "\n",
    "Шпаргалка по матричным производным:\n",
    "\n",
    "\n",
    "$$\\large \\begin{array}{rcl} \n",
    "\\frac{\\partial}{\\partial \\textbf{X}} \\textbf{X}^{\\text{T}} \\textbf{A} &=& \\textbf{A} \\\\\n",
    "\\frac{\\partial}{\\partial \\textbf{X}} \\textbf{X}^{\\text{T}} \\textbf{A} \\textbf{X} &=& \\left(\\textbf{A} + \\textbf{A}^{\\text{T}}\\right)\\textbf{X} \\\\\n",
    "\\frac{\\partial}{\\partial \\textbf{A}} \\textbf{X}^{\\text{T}} \\textbf{A} \\textbf{y} &=&  \\textbf{X}^{\\text{T}} \\textbf{y}\\\\\n",
    "\\frac{\\partial}{\\partial \\textbf{X}} \\textbf{A}^{-1} &=& -\\textbf{A}^{-1} \\frac{\\partial \\textbf{A}}{\\partial \\textbf{X}} \\textbf{A}^{-1} \n",
    "\\end{array}$$\n",
    "\n",
    "Продолжим:\n",
    "\n",
    "$$\\Large \\begin{array}{rcl} \\frac{\\partial \\mathcal{L}}{\\partial \\textbf{w}} &=& \\frac{\\partial}{\\partial \\textbf{w}} \\frac{1}{2n} \\left( \\textbf{y}^{\\text{T}} \\textbf{y} -2\\textbf{y}^{\\text{T}} \\textbf{X} \\textbf{w} + \\textbf{w}^{\\text{T}} \\textbf{X}^{\\text{T}} \\textbf{X} \\textbf{w}\\right) \\\\\n",
    "&=& \\frac{1}{2n} \\left(-2 \\textbf{X}^{\\text{T}} \\textbf{y} + 2\\textbf{X}^{\\text{T}} \\textbf{X} \\textbf{w}\\right)\n",
    "\\end{array}$$\n",
    "\n",
    "$$\\Large \\begin{array}{rcl} \\frac{\\partial \\mathcal{L}}{\\partial \\textbf{w}} = 0 &\\Leftrightarrow& \\frac{1}{2n} \\left(-2 \\textbf{X}^{\\text{T}} \\textbf{y} + 2\\textbf{X}^{\\text{T}} \\textbf{X} \\textbf{w}\\right) = 0 \\\\\n",
    "&\\Leftrightarrow& -\\textbf{X}^{\\text{T}} \\textbf{y} + \\textbf{X}^{\\text{T}} \\textbf{X} \\textbf{w} = 0 \\\\\n",
    "&\\Leftrightarrow& \\textbf{X}^{\\text{T}} \\textbf{X} \\textbf{w} = \\textbf{X}^{\\text{T}} \\textbf{y} \\\\\n",
    "&\\Leftrightarrow& \\textbf{w} = \\left(\\textbf{X}^{\\text{T}} \\textbf{X}\\right)^{-1} \\textbf{X}^{\\text{T}} \\textbf{y}\n",
    "\\end{array}$$\n",
    "\n",
    "Итак, имея в виду все определения и условия описанные выше, мы можем утверждать, опираясь на <a href=\"https://ru.wikipedia.org/wiki/%D0%A2%D0%B5%D0%BE%D1%80%D0%B5%D0%BC%D0%B0_%D0%93%D0%B0%D1%83%D1%81%D1%81%D0%B0_%E2%80%94_%D0%9C%D0%B0%D1%80%D0%BA%D0%BE%D0%B2%D0%B0\">теорему Маркова-Гаусса</a>, что оценка МНК является лучшей оценкой параметров модели, среди всех <i>линейных</i> и <i>несмещенных</i> оценок, то есть обладающей наименьшей дисперсией."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Метод максимального правдоподобия\n",
    "\n",
    "У читателя вполне резонно могли возникнуть вопросы: например, почему мы минимизируем среднеквадратичную ошибку, а не что-то другое. Ведь можно мнимизировать среднее абсолютное значение невязки или еще что-то. Единственное, что произойдёт в случае изменения минимизируемого значения, так это то, что мы выйдем из условий теоремы Маркова-Гаусса и наши оценки перестанут быть лучшими среди линейных и несмещенных. \n",
    "\n",
    "Давайте перед тем как продолжить, сделаем лирическое отступление, чтобы проиллюстрировать метод максимального правдоподобия на простом примере. \n",
    "\n",
    "Как-то после школы я заметил, что все помнят формулу этилового спирта. Тогда я решил провести эксперимент: помнят ли люди более простую формулу метилового спирта: $CH_3OH$. Мы опросили 400 человек и оказалось, что формулу помнят всего 117 человек. Разумно предположить, что вероятность того, что следующий опрошенный знает формулу метилового спирта – $\\frac{117}{400} \\approx 29\\%$. Покажем, что такая интуитивно понятная оценка не просто хороша, а еще и является оценкой максимального правдоподобия. \n",
    "\n",
    "Разберемся, откуда берется эта оценка, а для этого вспомним определение <a href=\"https://ru.wikipedia.org/wiki/%D0%A0%D0%B0%D1%81%D0%BF%D1%80%D0%B5%D0%B4%D0%B5%D0%BB%D0%B5%D0%BD%D0%B8%D0%B5_%D0%91%D0%B5%D1%80%D0%BD%D1%83%D0%BB%D0%BB%D0%B8\">распределения Бернулли</a>: случайная величина $X$ имеет распределение Бернулли, если она принимает всего два значения ($1$ и $0$ с вероятностями $\\theta$ и $1 - \\theta$ соответственно) и имеет следующую функцию распределения вероятности:\n",
    "\n",
    "$$\\large p\\left(\\theta, x\\right) = \\theta^x \\left(1 - \\theta\\right)^\\left(1 - x\\right),~x \\in \\left\\{0, 1\\right\\}$$\n",
    "\n",
    "Похоже, это распределение – то, что нам нужно, а параметр распределения $\\theta$ и есть та оценка вероятности того, что человек знает формулу метилового спирта. Мы проделали $400$ <i>независимых</i> экспериментов, обозначим их исходы как $\\textbf{x} = \\left(x_1, x_2, \\ldots, x_{400}\\right)$. Запишем *правдоподобие* наших данных (наблюдений), то есть вероятность наблюдать 117 реализаций случайной величины $X = 1$ и 283 реализации  $X = 0$:\n",
    "\n",
    "$$\\large p(\\textbf{x}; \\theta) = \\prod_{i=1}^{400} \\theta^{x_i} \\left(1 - \\theta\\right)^{\\left(1 - x_i\\right)} = \\theta^{117} \\left(1 - \\theta\\right)^{283}$$\n",
    "\n",
    "Далее будем максимизировать это выражение по $\\theta$, и чаще всего это делают не с правдоподобием $p(\\textbf{x}; \\theta)$, а с его логарифмом (применение монотонного преобразования не изменит решение, но упростит вычисления):\n",
    "\n",
    "$$\\large \\log p(\\textbf{x}; \\theta) = \\log \\prod_{i=1}^{400} \\theta^{x_i} \\left(1 - \\theta\\right)^{\\left(1 - x_i\\right)} = $$\n",
    "$$ \\large = \\log \\theta^{117} \\left(1 - \\theta\\right)^{283} =  117 \\log \\theta + 283 \\log \\left(1 - \\theta\\right)$$\n",
    "\n",
    "Теперь мы хотим найти такое значение $\\theta$, которое максимизирует правдоподобие, для этого мы возьмем производную по $\\theta$, приравняем к нулю и решим полученное уравнение:\n",
    "\n",
    "$$\\large  \\frac{\\partial p(\\textbf{x}; \\theta)}{\\partial \\theta} = \\frac{\\partial}{\\partial \\theta} \\left(117 \\log \\theta + 283 \\log \\left(1 - \\theta\\right)\\right) = \\frac{117}{\\theta} - \\frac{283}{1 - \\theta};$$\n",
    "\n",
    "\n",
    "$$\\large \\begin{array}{rcl} \n",
    "\\frac{117}{\\theta} - \\frac{283}{1 - \\theta} = 0 \\Rightarrow \\theta = \\frac{117}{400}\n",
    "\\end{array}.$$\n",
    "\n",
    "Получается, что наша интуитивная оценка – это и есть оценка максимального правдоподобия.  Применим теперь те же рассуждения для задачи линейной регрессии и попробуем выяснить, что лежит за среднеквадратичной ошибкой. Для этого нам придется посмотреть на линейную регрессию с вероятностной точки зрения. Модель, естественно, остается такой же:\n",
    "\n",
    "$$\\large \\textbf y = \\textbf X \\textbf w + \\epsilon,$$\n",
    "\n",
    "но будем теперь считать, что случайные ошибки берутся из центрированного <a href=\"https://ru.wikipedia.org/wiki/%D0%9D%D0%BE%D1%80%D0%BC%D0%B0%D0%BB%D1%8C%D0%BD%D0%BE%D0%B5_%D1%80%D0%B0%D1%81%D0%BF%D1%80%D0%B5%D0%B4%D0%B5%D0%BB%D0%B5%D0%BD%D0%B8%D0%B5\">нормального распределения</a>:\n",
    "\n",
    "$$\\large \\epsilon_i \\sim \\mathcal{N}\\left(0, \\sigma^2\\right)$$\n",
    "\n",
    "Перепишем модель в новом свете:\n",
    "\n",
    "$$\\large \\begin{array}{rcl} \n",
    "y_i &=& \\sum_{j=0}^m w_j X_{ij} + \\epsilon_i \\\\\n",
    "&\\sim& \\sum_{j=0}^m w_j  X_{ij} + \\mathcal{N}\\left(0, \\sigma^2\\right) \\\\\n",
    "p\\left(\\textbf{y} \\mid \\textbf X; \\textbf{w}\\right) &=& \\mathcal{N}\\left(\\sum_{j=0}^m w_j X_{ij}, \\sigma^2\\right)\n",
    "\\end{array}$$\n",
    "\n",
    "Так как примеры берутся независимо (ошибки не скоррелированы – одно из условий теоремы Маркова-Гаусса), то полное правдоподобие данных будет выглядеть как произведение функций плотности $p\\left(y_i\\right)$. Рассмотрим логарифм правдоподобия, что позволит нам перейти от произведения к сумме:\n",
    "\n",
    "$$\\Large \\begin{array}{rcl} \n",
    "\\log p\\left(\\textbf{y} \\mid \\textbf X; \\textbf{w}\\right) &=& \\log \\prod_{i=1}^n \\mathcal{N}\\left(\\sum_{j=0}^m w_j X_{ij}, \\sigma^2\\right) \\\\\n",
    "&=& \\sum_{i=1}^n \\log \\mathcal{N}\\left(\\sum_{j=0}^m w_j X_{ij}, \\sigma^2\\right) \\\\\n",
    "&=& -\\frac{n}{2}\\log 2\\pi\\sigma^2 -\\frac{1}{2\\sigma^2} \\sum_{i=1}^n \\left(y_i - \\textbf{w}^{\\text{T}} \\textbf{x}_i\\right)^2\n",
    "\\end{array}$$\n",
    "\n",
    "Мы хотим найти гипотезу максимального правдоподобия, т.е. нам нужно максимизировать выражение $p\\left(\\textbf{y} \\mid \\textbf X; \\textbf{w}\\right)$ по вектору $\\textbf w$, получив при этом $\\textbf{w}_{\\text{ML}}$. Обратите внимание, что при максимизации функции по какому-то параметру можно выкинуть все члены, не зависящие от этого параметра:\n",
    "\n",
    "$$\\Large \\begin{array}{rcl} \n",
    "\\textbf{w}_{\\text{ML}} &=& \\arg \\max_{\\textbf w} p\\left(\\textbf{y}\\mid \\textbf X; \\textbf{w}\\right) = \\arg \\max_{\\textbf w} \\log p\\left(\\textbf{y}\\mid \\textbf X; \\textbf{w}\\right)\\\\\n",
    "&=& \\arg \\max_{\\textbf w} -\\frac{n}{2}\\log 2\\pi\\sigma^2 -\\frac{1}{2\\sigma^2} \\sum_{i=1}^n \\left(y_i - \\textbf{w}^{\\text{T}} \\textbf{x}_i\\right)^2 \\\\\n",
    "&=& \\arg \\max_{\\textbf w} -\\frac{1}{2\\sigma^2} \\sum_{i=1}^n \\left(y_i - \\textbf{w}^{\\text{T}} \\textbf{x}_i\\right)^2 \\\\\n",
    "&=&  \\arg \\max_{\\textbf w} \\mathcal{L}\\left(\\textbf X, \\textbf{y}, \\textbf{w} \\right)\n",
    "\\end{array}$$\n",
    "\n",
    "Таким образом, мы увидели, что максимизация правдоподобия данных – это то же самое, что и минимизация среднеквадратичной ошибки (при справедливости указанных выше предположений). Получается, что именно такая функция стоимости является следствием того, что ошибка распределена нормально, а не как-то по-другому."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Разложение ошибки на смещение и разброс (Bias-variance decomposition)\n",
    "\n",
    "Поговорим немного о свойствах оценки, полученной линейной регрессией. В свете предыдущего пункта мы выяснили, что:\n",
    "- истинное значение целевой переменной складывается из некоторой детерминированной функции $f(\\textbf{x})$ и случайной ошибки $\\epsilon$: $y = f\\left(\\textbf{x}\\right) + \\epsilon$;\n",
    "- ошибка распределена нормально с центром в нуле и некоторым разбросом: $\\epsilon \\sim \\mathcal{N}\\left(0, \\sigma^2\\right)$;\n",
    "- истинное значение целевой переменной тоже распределено нормально: $y \\sim \\mathcal{N}\\left(f\\left(\\textbf{x}\\right), \\sigma^2\\right)$\n",
    "- мы пытаемся приблизить детерминированную, но неизвестную функцию $f\\left(\\textbf{x}\\right)$ линейной функцией от регрессоров $\\widehat{f}\\left(\\textbf{x}\\right)$, которая, в свою очередь, является точечной оценкой функции $f$ в пространстве функций (точнее, мы ограничили пространство функций параметрическим семейством линейных функций), т.е. случайной переменной, у которой есть среднее значение и дисперсия.\n",
    "\n",
    "Тогда ошибка в точке $\\textbf{x}$ раскладывается следующим образом:\n",
    "\n",
    "$$\\large \\begin{array}{rcl} \n",
    "\\text{Err}\\left(\\textbf{x}\\right) &=& \\mathbb{E}\\left[\\left(y - \\widehat{f}\\left(\\textbf{x}\\right)\\right)^2\\right] \\\\\n",
    "&=& \\mathbb{E}\\left[y^2\\right] + \\mathbb{E}\\left[\\left(\\widehat{f}\\left(\\textbf{x}\\right)\\right)^2\\right] - 2\\mathbb{E}\\left[y\\widehat{f}\\left(\\textbf{x}\\right)\\right] \\\\\n",
    "&=& \\mathbb{E}\\left[y^2\\right] + \\mathbb{E}\\left[\\widehat{f}^2\\right] - 2\\mathbb{E}\\left[y\\widehat{f}\\right] \\\\\n",
    "\\end{array}$$\n",
    "\n",
    "Для наглядности опустим обозначение аргумента функций. Рассмотрим каждый член в отдельности, первые два расписываются легко по формуле $\\text{Var}\\left(z\\right) = \\mathbb{E}\\left[z^2\\right] - \\mathbb{E}\\left[z\\right]^2$:\n",
    "\n",
    "$$\\large \\begin{array}{rcl} \n",
    "\\mathbb{E}\\left[y^2\\right] &=& \\text{Var}\\left(y\\right) + \\mathbb{E}\\left[y\\right]^2 = \\sigma^2 + f^2\\\\\n",
    "\\mathbb{E}\\left[\\widehat{f}^2\\right] &=& \\text{Var}\\left(\\widehat{f}\\right) + \\mathbb{E}\\left[\\widehat{f}\\right]^2 \\\\\n",
    "\\end{array}$$\n",
    "\n",
    "Пояснения:\n",
    "\n",
    "$$\\large \\begin{array}{rcl} \n",
    "\\text{Var}\\left(y\\right) &=& \\mathbb{E}\\left[\\left(y - \\mathbb{E}\\left[y\\right]\\right)^2\\right] \\\\\n",
    "&=& \\mathbb{E}\\left[\\left(y - f\\right)^2\\right] \\\\\n",
    "&=& \\mathbb{E}\\left[\\left(f + \\epsilon - f\\right)^2\\right] \\\\\n",
    "&=& \\mathbb{E}\\left[\\epsilon^2\\right] = \\sigma^2\n",
    "\\end{array}$$\n",
    "\n",
    "$$\\large \\mathbb{E}[y] = \\mathbb{E}[f + \\epsilon] = \\mathbb{E}[f] + \\mathbb{E}[\\epsilon] = f$$\n",
    "\n",
    "И теперь последний член суммы. Мы помним, что ошибка и целевая переменная независимы друг от друга:\n",
    "\n",
    "$$\\large \\begin{array}{rcl} \n",
    "\\mathbb{E}\\left[y\\widehat{f}\\right] &=& \\mathbb{E}\\left[\\left(f + \\epsilon\\right)\\widehat{f}\\right] \\\\\n",
    "&=& \\mathbb{E}\\left[f\\widehat{f}\\right] + \\mathbb{E}\\left[\\epsilon\\widehat{f}\\right] \\\\\n",
    "&=& f\\mathbb{E}\\left[\\widehat{f}\\right] + \\mathbb{E}\\left[\\epsilon\\right] \\mathbb{E}\\left[\\widehat{f}\\right]  = f\\mathbb{E}\\left[\\widehat{f}\\right]\n",
    "\\end{array}$$\n",
    "\n",
    "Наконец, собираем все вместе:\n",
    "\n",
    "$$\\large \\begin{array}{rcl} \n",
    "\\text{Err}\\left(\\textbf{x}\\right) &=& \\mathbb{E}\\left[\\left(y - \\widehat{f}\\left(\\textbf{x}\\right)\\right)^2\\right] \\\\\n",
    "&=& \\sigma^2 + f^2 + \\text{Var}\\left(\\widehat{f}\\right) + \\mathbb{E}\\left[\\widehat{f}\\right]^2 - 2f\\mathbb{E}\\left[\\widehat{f}\\right] \\\\\n",
    "&=& \\left(f - \\mathbb{E}\\left[\\widehat{f}\\right]\\right)^2 + \\text{Var}\\left(\\widehat{f}\\right) + \\sigma^2 \\\\\n",
    "&=& \\text{Bias}\\left(\\widehat{f}\\right)^2 + \\text{Var}\\left(\\widehat{f}\\right) + \\sigma^2\n",
    "\\end{array}$$\n",
    "\n",
    "Итак, мы достигли цели всех вычислений, описанных выше, последняя формула говорит нам, что ошибка прогноза любой модели вида $y = f\\left(\\textbf{x}\\right) + \\epsilon$ складывается из:\n",
    "\n",
    "- квадрата смещения: $\\text{Bias}\\left(\\widehat{f}\\right)$ – средняя ошибка по всевозможным наборам данных;\n",
    "- дисперсии: $\\text{Var}\\left(\\widehat{f}\\right)$ – вариативность ошибки, то, на сколько ошибка будет отличаться, если обучать модель на разных наборах данных;\n",
    "- неустранимой ошибки: $\\sigma^2$.\n",
    "\n",
    "Если с последней мы ничего сделать не можем, то на первые два слагаемых мы можем как-то влиять. В идеале, конечно же, хотелось бы свести на нет оба этих слагаемых (левый верхний квадрат рисунка), но на практике часто приходится балансировать между смещенными и нестабильными оценками (высокая дисперсия).\n",
    "\n",
    "<img src=\"../../img/bvtf.png\" width=\"480\">\n",
    "\n",
    "Как правило, при увеличении сложности модели (например, при увеличении количества свободных параметров) уменьшается дисперсия (разброс) оценки, но оценка становится смещенной. Из-за того что тренировочный набор данных полностью запоминается вместо обобщения, небольшие изменения приводят к неожиданным результатам (переобучение). Если же модель слабая, то она не в состоянии выучить закономерность, в результате выучивается что-то другое, смещенное относительно правильного решения.\n",
    "\n",
    "<img src=\"../../img/biasvariance.png\" width=\"480\">\n",
    "\n",
    "Теорема Маркова-Гаусса как раз утверждает, что МНК-оценка параметров линейной модели является самой лучшей в классе несмещенных линейных оценок, то есть с наименьшей дисперсией. Это значит, что если существует какая-либо другая несмещенная модель $g$ тоже из класса линейных моделей, то мы можем быть уверены, что $Var\\left(\\widehat{f}\\right) \\leq Var\\left(g\\right)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Регуляризация линейной регрессии\n",
    "\n",
    "Иногда бывают ситуации, когда мы намеренно увеличиваем смещенность модели ради ее стабильности, т.е. ради уменьшения дисперсии модели $\\text{Var}\\left(\\widehat{f}\\right)$. Одним из условий теоремы Маркова-Гаусса является полный столбцовый ранг матрицы $\\textbf{X}$. В противном случае решение МНК $\\textbf{w} = \\left(\\textbf{X}^{\\text{T}} \\textbf{X}\\right)^{-1} \\textbf{X}^{\\text{T}} \\textbf{y}$ не существует, т.к. не будет существовать обратная матрица $\\left(\\textbf{X}^{\\text{T}} \\textbf{X}\\right)^{-1}.$ Другими словами, матрица $\\textbf{X}^{\\text{T}} \\textbf{X}$ будет вырождена. Такая задача называется <a href=\"https://ru.wikipedia.org/wiki/%D0%9A%D0%BE%D1%80%D1%80%D0%B5%D0%BA%D1%82%D0%BD%D0%BE_%D0%BF%D0%BE%D1%81%D1%82%D0%B0%D0%B2%D0%BB%D0%B5%D0%BD%D0%BD%D0%B0%D1%8F_%D0%B7%D0%B0%D0%B4%D0%B0%D1%87%D0%B0\">некорректно поставленной</a>. Задачу нужно скорректировать, а именно, сделать матрицу $\\textbf{X}^{\\text{T}}\\textbf{X}$ невырожденной, или регулярной (именно поэтому этот процесс называется регуляризацией). Чаще в данных мы можем наблюдать так называемую <i>мультиколлинеарность</i> — когда два или несколько признаков сильно коррелированы, в матрице $\\textbf{X}$ это проявляется в виде \"почти\" линейной зависимости столбцов. Например, в задаче прогнозирования цены квартиры по ее параметрам \"почти\" линейная зависимость будет у признаков \"площадь с учетом балкона\" и \"площадь без учета балкона\". Формально для таких данных матрица $\\textbf{X}^{\\text{T}} \\textbf{X}$ будет обратима, но из-за мультиколлинеарности у матрицы $\\textbf{X}^{\\text{T}} \\textbf{X}$ некоторые собственные значения будут близки к нулю, а в обратной матрице $\\left(\\textbf{X}^{\\text{T}} \\textbf{X}\\right)^{-1}$ появятся экстремально большие собственные значения, т.к. собственные значения обратной матрицы – это $\\frac{1}{\\lambda_i}$. Итогом такого шатания собственных значений станет нестабильная оценка параметров модели, т.е. добавление нового наблюдения в набор тренировочных данных приведёт к совершенно другому решению. Иллюстрации роста коэффициентов вы найдете в <a href=\"https://habrahabr.ru/company/ods/blog/322076/\">одном из наших прошлых постов</a>. Одним из способов регуляризации является <a href=\"https://ru.wikipedia.org/wiki/%D0%9C%D0%B5%D1%82%D0%BE%D0%B4_%D1%80%D0%B5%D0%B3%D1%83%D0%BB%D1%8F%D1%80%D0%B8%D0%B7%D0%B0%D1%86%D0%B8%D0%B8_%D0%A2%D0%B8%D1%85%D0%BE%D0%BD%D0%BE%D0%B2%D0%B0\">регуляризация Тихонова</a> (a.k.a. L2-регуляризация), которая в общем виде выглядит как добавление нового члена к среднеквадратичной ошибке: \n",
    "\n",
    "$$\\Large \\begin{array}{rcl} \n",
    "\\mathcal{L}\\left(\\textbf{X}, \\textbf{y}, \\textbf{w} \\right) &=& \\frac{1}{2n} \\left\\| \\textbf{y} - \\textbf{X} \\textbf{w} \\right\\|_2^2 + \\left\\|\\Gamma \\textbf{w}\\right\\|^2\\\\\n",
    "\\end{array}$$\n",
    "\n",
    "Часто матрица Тихонова выражается как произведение некоторого числа на единичную матрицу: $\\Gamma = \\frac{\\lambda}{2}  \\textbf{E}$. В этом случае задача минимизации среднеквадратичной ошибки становится задачей с ограничением на $L_2$-норму. Если продифференцировать новую функцию стоимости по параметрам модели, приравнять полученную функцию к нулю и выразить $\\textbf{w}$, то мы получим точное решение задачи.\n",
    "\n",
    "$$\\Large \\begin{array}{rcl} \n",
    "\\textbf{w} &=& \\left(\\textbf{X}^{\\text{T}} \\textbf{X} + \\lambda \\textbf{E}\\right)^{-1} \\textbf{X}^{\\text{T}} \\textbf{y}\n",
    "\\end{array}$$\n",
    "\n",
    "Такая регрессия называется гребневой регрессией (ridge regression). А гребнем является как раз диагональная матрица, которую мы прибавляем к матрице $\\textbf{X}^{\\text{T}} \\textbf{X}$, в результате получается гарантированно регулярная матрица.\n",
    "\n",
    "<img src=\"../../img/ridge.png\">\n",
    "\n",
    "Такое решение уменьшает дисперсию, но становится смещенным, т.к. минимизируется также и норма вектора параметров, что заставляет решение сдвигаться в сторону нуля. На рисунке ниже на пересечении белых пунктирных линий находится МНК-решение. Голубыми точками обозначены различные решения гребневой регрессии. Видно, что при увеличении параметра регуляризации $\\lambda$ решение сдвигается в сторону нуля.\n",
    "\n",
    "<img src=\"../../img/l2.png\">\n",
    "\n",
    "Советуем обратиться в <a href=\"https://habrahabr.ru/company/ods/blog/322076/\">наш прошлый пост</a> за примером того, как $L_2$-регуляризация справляется с проблемой мультиколлинеарности, а также чтобы освежить в памяти еще несколько интерпретаций регуляризации."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
